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Although  it  is  generally  recognized  that  texture  i mares 
cor. uain  statistical,  spectral  and  structural  domain  information, 
the  use  of  spectral  information  alone  can  he  quite  effective  in 
the  texture-image  analysis  studies  such  as  texture  discrimination 
and  segmentation.  Bajcsy  and  Liberman  [i]  expressed  the  power 
spectrum  in  polar  coordinates,  then  integrate  over  r  and  IQ,  to 
obtain  the  two  one-dimensional  functions.  The  location  cf  peaks 
in  these  functions  indicates  prominant  texture  coarseness  and 
directionality.  Wes  ska  et.  ai.  [  2]  integrated  the  power  spectrum 
within  16  spatial  frequency  zones  which  were  combinations  of  four 
1-octave  frequency  ranges  and  four  45°  orientation  sectors.  They 
also  computed  sight  "contrast"  measures  based  on  the  cooccurrence 
matrix,  and  obtained  better  discrimination  than  with  the  power 
spectrum  measures.  Laws  £3]  computed  a  number  of  energy  aeasurt  :■ 
by  filtering  the  texture  with  sets  of  small  linear  operators, 
then  squaring  and  summing  the  output  of  each  filter.  He  report'  : 
better  discrimination  with  the  energy  than  with  the  cooccurrence 
measures . 

A  fundamental  problem  with  the  power  spectrum  analysis  is  the 
computational  accuracy  and  computational  complexity.  For  texture 
study,  accurate  power  spectrum  must  be  computed  from  the  small 
image  segments.  In  this  case,  the  two-dimensional  Fourier  analyst 
cannot  provide  sufficient  accuracy  as  the  Fourier  analysis  is  mere- 
accurate  with  a  large  number  of  pixels.  The  two-dimensional 
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maximum  entropy  spectral  analysis,  however,  is  very  suitable  for 
a  small  number  of  pixels.  The  computational  complexity  has  been 
a  drawback  in  using  the  two-dimensional  maximum  entropy  spectral 
estimation  procedures.  Recently,  Lim  and  Wlalik  [^4,  5»  have 
proposed  an  efficient  iterative  algorithm  for  the  two -dimensions] 
maximum  entropy  power  spectrum  estimation  suitable  for  the  minic.n:: 
puter  implementation.  Their  method  is  adapted  and  generalized  f  r 
use  in  our  PDF  11/45  minicomputer  for  the  texture-image  analysis. 

A  three-dimensional  graphics  software  is  developed  for  the  spec¬ 
tral  display  at  the  different  viewing  positions.  Extensive  com¬ 
puter  results  on  the  spectral  analysis  of  texture  images  are  also 
reported . 

II.  Two-Dimensional  Power  Spectrum  Estimation 

To  obtain  the  power  spectrum  of  a  two-dimensional  signal, 
the  direct  method  is  to  calculate  the  two-dimensional  Fourier 


transform  of  the  autocorrelation  function,  i.e.s 

r>»  op 

Fv(w 


w. 


)  =  Z  Z 

The  following  notations  will  be  used  in  the  report: 


(1 


x(n^,  n£ ) :  A  2-D  random  signal  whose  power  spectrum  we  wish 
to  estimate. 


Vv 

n2h 

Autocorrelation  function  of  x(n^,  n2 ) 

Rx(nl, 

n2): 

An  estimate  of  R  (r.^,  ) 

Vv 

w2)s 

Power  spectrum  of  xCn^,  ) 

A  (n1 , 

n2): 

Autocorrelation  function  whose  power  spectrum 
l/?x(wi»  w2) 

A 

: 

A  set  of  points  (n^,  )  for  which  Rx(n^,  n2 ) 

known. 

F 

: 

Discrete  time  Fourier  Transform 

F_1 

. 

Inverse  discrete  time  Fourier  Transform 
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From  (1),  it  is  important  to  note  that  the  determination  o: 
the  power  spectral  density  entails  complete  knowledge  of  the 
generally  infinite  extent  autocorrelation  function.  For  the 
finite  signal,  this  method  is  proved  to  have  poor  resolution  due 
to  the  truncated  and  s-^mpled  autocorrelation  function  set. 


There  are  various  techniques  to  estimate  the  power  spectruri 
for  the  one-dimensional  signal  M-  One  technique  that  has  beer, 
recognized  as  the  best  due  to  its  high  resolution  is  the  Maximum 
Entropy  Method  (MEM).  The  basic  idea  of  this  approach  is  to  extr 
polate  the  autocorrelation  function  of  a  random  process  by  max  it.:.! 


the  entropy  H  of  the  corresponding  probability  density  function: 

H  =  j  log  Fx(w)  dw  (2) 

UVr-TT 

where  P  (w)  is  the  power  spectrum  density.  The  characteristics  ■„ 
the  maximum  entropy  method  are  equivalent  to  the  autoregressive  ' 
signal  modeling  £sj  which  requires  solving  a  set  of  linear  equa'- i 
for  the  filter  coefficients.  This  can  be  expressed  as: 

A  / 

/  \  t  (  r> 

?  (w)  = 


M  -;«0  2 

1  +  H  ak  e  J 


and  the  filter  coefficients  are  obtained  by  solving  the  norma. 


equations,  R  (i) 


-  21  a k  R(1-k) 


There  are  different  algorithms  proposed  to  find  the  solutions  of 


the  normal  equations.  The  most  efficient  ones  are  the  Levinson 
recursive  algorithm  and  the  Burg  recursive  algorithm  [ ?]  •  But 
for  the  two-dimensional  case,  the  problem  is  different  since  the 
normal  equations  become  a  highly  nonlinear  problem  [4].  For  the 


general  form  of  the  two-dimensional  case, 


(5) 


2  Iai.i  V1"1'  “*i)  *  Rx  <r,  =) 


Ci,j)fB  for  ( r>  s)  ^  E 

Mere  the  set  B  consists  of  all  points  where  the  filter  mask  has 

non-zero  values,  and  the  power  spectrum  obtained  from  a. .  is  given 
^  ^  ^ 

by  Px<Wl'  "2)  =  l l6> 

In  this  case,  we  can  see  from  (5)  that  the  size  of  independent 
values  of  R  (n, ,  n9 )  requirea  to  solve  the  above  sex  of  equations  is 
FL  eater  than  the  size  of  the  filter  mask.  For  example,  Fig.  1(a) 
?rnv,\:  the  autoregressive  filter  mask  size  as  3x2 ,  ana  Fig.  1(c) 
shown  a  larger  size  of  independent  values  <-f  K  (rM  ,  r... )  reauired  to 
solve  for  a-  .  in  rig.  1(a)  by  equation  (5;. 


Fig.  1 

Clearly,  the  number  of  correlation  points  needed  is  greater  than 
the  number  of  filter  coefficients.  Since  the  estimated  power 
spectrum  given  by  (6)  is  completely  determined  by  the  coefficients 


alone,  it  does  not  posses  enough  degrees  of  freedom  to  solve  for 
the  spectrum.  Therefore,  the  normal  equations  are  not  linear  as 
in  the  one-dimensional  case.  Many  methods  have  been  proposed  to 
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t'x  tend  the  Levinson's  and  the  Burg's  algorithms  in  two  dimension, 
however,  those  algorithms  are  not  computationally  attractive,  and 
: hero  is  no  guarantee  that  a  solution  or  an  approximate  solution  can 
re  retained.  For  instance,  Burg  j^l 3 J  has  proposed  an  iterative 
solution  which  requires  the  inversion  of  a  matrix  in  each  iteration 
•where  the  dimension  of  the  matrix  is  of  the  order  of  the  number  of 
the  riven  autocorrelation  points.  No  experimental  results  using 
:iiis  technique  have  yet  been  reported.  Wernecke  and  D'Addario  [lhj 
have  proposed  a  scheme  in  which  an  attempt  is  made  to  numerically 
maximize  the  entropy.  The  maximization  is  done  by  continuously 
adjusting  the  power  spectrum  estimate  and  evaluating  the  expressions 
t'or  the  entropy  and  its  gradient.  The  procedure  is  computationally 
expensive  and  is  not  guaranteed  to  have  a  solution.  Woods  [ll] 
expresses  the  Maximum  Entropy  Method  as  a  power  series  in  the 
frequency  domain  and  attempts  to  approximate  the  ME  FS  estimate  by 
truncating  the  power  series  expansion.  Even  though  such  an  approach 
has  some  computational  advantages  relative  to  others,  the  method  is 
restricted  to  the  class  of  signals  for  which  the  power  series  expansi  >n 
is  possible. 

Based  on  the  reason  that  the  closed  form  solution  of  the  two- 
dimensional  ME  method  is  hard  to  obtain,  Lim  and  Malik  developed  a 
now  iterative  algorithm,  using  adaptive  filtering  concept.  This 
u 1 rorithm  can  correctly  estimate  the  true  spectrum  and  is  compute - 
’  :  i roily  simple  due  to  the  utilization  of  Fast  Fourier  Transform 
( fTT '• .  The  basic  idea  of  this  algorithm  is  on  the  notion  that  the 
r’  ''•■n  correlation  points  in  region  A  is  consistent  ana  the  con\i> 


-  o  - 


ponding  coefficient  shoulc  be  zero  outside  region  A,  and  proceed 
this  iteration  repeatedly  until  the  optimal  solution  is  obtained. 

A 

Thai  is,  given  RJ£(n1,  N?)  for  (m  ,  n2)£  A, determine  P  (v^,  w2  ] 

A 

such  that  Px(wi_»  v;0)  has  the  form 


Px  (V  V 


1 


(7) 


and 


Rx ( n i *  n2)  =  F_1  [px(wlf  w2)J  for  (n1 ,  n£)  £  A 
A  sin-pie  flowchart  is  shown  in  Fig.  2.  We  begin  with  some  initial 
estimate  of  A(n^,  n2 ) ,  obtain  the  corresponding  correlation  function, 
correct  the  resulting  correlation  function  for  (n^,  n2 )  £  A  with  the 
Known  Rx(n1 ,  n2 ) ,  obtain  the  corresponding  A(n^,  n2 )  from  the 
correct  correlation  function,  and  then  replace  the  resulting  A(n.,n  ) 
"Uh  o  for  (n  ,  n2)  A.  This  completes  one  iteration  and  the  correct' a 
A  (n^,  n2)  is  a  new  estimate  of  A (n1 ,  n2). 

initial  estimate  of  (n^,  n2) 

I 


->•  R  (n1(  n2)  =  F 


-1 


\ 


Ffxr, 


•VM]  ] 


correct  R  (n^,  n2)  with  Rx(^1,  n2)  for  (n1,  n?)  £  A 


I 


L 


A  (ni’  n2}  =  F  *[  FfVni-n^]j 

Y 

?\  (n1,  n2)  =  o  for  (n1>  n?)  ^  A 

I 

Px(wl’  W2}  =  F  [Ry  (nl>  n2}] 


Fig.  2 


However ,  to  prevent  the  zero  crossing  problem,  when  taking  xhe  inverse 
01  P  ’  n2^J  an<^  ^  [  ^y^ni*  ng)J  due  "t0  ”^-r;e  correction  of 

1'.,  (n.j  ,  n2 )  and  7\  (n^,  ) ,  and  also  to  keep  the  correlation  function 

to  be  positive  definite,  they  modify  the  procedure  by  linearly  inter¬ 
polating  some  parameters  to  prevent  divergence  of  the  error  and  to 
increase  the  rate  of  convergence  [h]  .  Thus  they  have  added  some 
constraints  to  make  the  practical  algorithm  as  shown  in  Fig.  3  . 

III.  Examples 

The  following  are  some  examples  of  using  Lim  and  Malik's 
algorithm  which  is  implemented  in  our  PDF  11/45  minicomputer.  The 
input;  signals  are  two-dimensional  autocorrelation  function  originated 
from  sinusoids  buried  in  white  noise,  so  that  the  correlation  func¬ 
tion  given  has  the  form  of 

M  _ 

Rx(n1  ,  n2)  =  O’*  5  («■]_,  n2)  +  21  ai“  costv;^  n.±+  wi2  r.,  )  (c) 

1-I 

2 

where  <T  is  the  white  noise  power,  }.]  is  the  number  of  sinusoids, 
a  j  ‘  /Z  is  the  power  of  the  ith  sinusoid,  and  and  represent  the 

frequency  of  the  ith  sinusoid.  Fig.  4  -  Fig.  6  are  the  cases  for 
1,  2,  3  sinusoids  respectively,  which  correspond  to  the  input  data 
snowii  in  Tables  1,  2,  3-  From  these  results,  we  can  see  this  iterativ 
procedure  can  easily  predict  the  true  spectrum  although  the  matrix 
of  autocorrelation  function  is  not  very  large. 

However,  in  Lim  and  Malik's  examples  and  the  examples  shown  in 
Pigs.  4-6,  the  expression  used  in  the  autocorrelation-function 

'  filiation  is  Eq.  (8)  which  is  an  analytical  form.  This  expression 
in  based  on  the  consideration  of  conxinuous  and  infinite  sinusoids. 

So  its  values  are  completely  symmetric.  For  the  practical  two- 
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GIVEN  R.{n,.nt).w(ni.n«).  DPT  LENCTH  N'.  k  =  O  S 
a.  =  0 ,f.=  0.  c.  =  IO-* 

INITIAL  ESTIMATE: 

Ry(ni,nt)  =  R.<0.0)  a(n,.n,) 

A*(n,.nt)  =  <("!■"») 


Fig.  3  A  detailed  flowchart  of  the  Lim- 
Malik  iterative  algorithm  for  2-d 
ME  PSE  implemented  in  the  report. 
UEEE  ASSP  Trans.  June  1981) 


•  •  i i ; ‘  ; :s lonui  signal  ,  xz.is  is  no  longer  tne  case.  le  general  izo  z;,i.s 

•  !.<:•'•  sithai  so  that  ix  can  ue  titteu  for  any  two-dimensional  signal , 


;he  well  known 


irmula  lor  the  calculation 


^nage  caxa. 


.c  .’”\ri.0 


I  N-n, 

w  £  £  xik’1  ;  xU™i'  (  "V 

La  expression  can  fix  the  real  autocorrelation 


:  -..no  •  ion,  we  have  used  ix  xo  generate  a  simulated  autocorrelation 

■  ion  of  the  two-dimensional  sinusoids,  and  the  spectrum  estimated 
•  •'  a  most  the  same  as  that  calculated  from  kq.  (£}. 


section,  we  will  ore sent  the  results  of 


two-dimens: 


a  x  i  mu: 


ropy  y.ethoa  lower  Spectrum 


:i:r.aticn  oi  tne  real  texture 


data  are  taken  from  the  U.  S. 


cata  case.  Ine  size 


■  p-ich  image  is  o-rx69.  she  original  textures  are  shown  in  Fig.  ?. 
i  pictures  will  reappear  but  5  times  larger  in  Fig.  6a  -  Fig.  16a 

with  its  corresponding  estimated  spectrum  shown  in  Fig.  6b  -  Fig.  lob 
rig.  6c  -  Fig.  16c  and  Fig.  cd  -  Fig.  l6d. 

Fig.  6b  shows  that  the  picture  has  one  main  frequency  component 
tv'ar  (0.01,  0.22),  one  small  dc  component  and  one  frequency  at 
i.  0 .  ,  0.5)  and  other  very  small  ripples.  Because  of  the  three- 
dimensional  display,  Fig.  be  and  Fig.  fid  cannot  accurately  indicate 
rhfi  (0.01,  0.22)  point,  but  we  can  see  the  entire  power  spectrum 
Cure , ion  distribution  in  Fig.  6c  and  Fig.  fid.  Fig.  96  shows  that 
ti;c  main  frequency  component  is  approximately  equal  xo  (0.02,0.126) 
and  ther  small  ripples.  Comparing  Figs.  8  ana  9,  we  can  see  clearly 
i  difference  in  power  spectrum  for  different  textures. 

Fig.  10  shows  that  the  texture  contains  a  main  frequency 
'Tient  around  (0.0,  0.0)  and  other  small  amplitude  components. 


J 


wen  the  corresponding  picture  (upper-  ri gh c  Sc; 

•n.  wee  that  it  contains  mostly  high  grey  level  :..i 
i  it  has  a  high  ac  compcr.ont.  By  taking  a  lo;.: 

•  i.-i  i  rodict  the  other  small  amplitude  irea uer.ov  t . 


;  n  0  Y/i:  in 


-■ •  1  j  union  is  a  r ec or. s tr u c 


—  *  11  .  o.  -  1 ..  .  *  *.  o ...  w  _l  ^  -j  o  f  . .  e  l.  ij.  : .  e  j .  .  _ . 


.ted.  is  tr.e  c  on. c  s  1 1  ion 


i  .  lid  and  I- is. 


:i,:»  star  \<j.  Oo,  ^ .  u;. ,  is  anecteu  'ey  the  contribution  of  rig.  12. 


■■•ten  we 


au  tiOCoxTc.it  oiC:.  x.  one  ,.oo  n.a  .inn  m  longer 


I  mgth,  tms  term  can  ce  recov; 


oiearr:  demonstrates 


o  .  o  onencm.encn. 


prove  the  accuracy  c;  the  two -dimensional  '.'.1 


-  I  arn- ithm,  three  sets  of  "nearly  periodic"  texture  data  were  tested. 

-  rr  in  Figs.  14-16,  tne  results  obtained  are  satisfactory  and  appear 
‘ v  indicate  truly  the  real  rower  spectrum. 

TV.  discussion 

For  a  real  two-dimensional  signal,  this  algorithm  can  accurately 
(relict  the  main  frequency  components  with  a  moderate  size  autccorre- 
j alien  function  (ACF)  matrix  and  short  DFT  length.  However,  for  the 
"'her-  frequency  terms  with  smaller  amplitude,  it  must  take  a  larger 
Adi-  matrix  and  longer  DFT  length  to  discriminate,  since  the  larger 
Adi'  matrix  will  add  more  information  about  the  signal.  But  this 
will  cake  too  much  computational  time.  Also  the  spectral  peaks  that 


1  cry  cicpfc  cannot  ce  resolvec 


:his  algorithm  with  the  moderate 


i/H  ACF  matrix  and  DFT  length.  To  understand  and  to  improve  the 
p<jr  !  ral  resolution  capability  of  the  algorithm,  we  will  analyze  the 
.  We  are  now  continuing  the  research  by  generating  an 


It  ue  spectrum 
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Fig.  7  The  original  test  texture  data 
(taken  fran  USC  data  base  ) .  Each  data 
format  is  64x64.  The  right  one  in 
second  row  is  reconstructed  from  the  left 
and  center  pictures  of  the  second  row. 


rig.  8(a)  The  original  test  data.  Fig.  8(b)  The  contour  map  with  JdZ3=3  with 

main  frequency  around  (0.01,  0.22). 


Fig.  8(d)  Three-dimensional  display  at 
(0.5, 0.5)  point  of  view. 


Fig.  8(c)  Three-dimensional  display  at 
(0.0, 0.0)  point  of  view. 
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Fig.  9(a)  The  original  test  data.  Fig.  9(b)  The  contour  map  with  <ldB=3  with 

main  frequency  around  (0.02,  0.126). 


Fig.  9(c)  Three-dimensional  display  at  Fig.  9(d)  Three-dimensional  display  au 

(0.0, 0.0)  point  of  view.  (0.5, 0.5)  point  of  view. 
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Fin.  lu(a)  The  original  test  data. 


Fig.  10(b)  The  contour  map  with  <4d3=3 
vath  main  frequency  around  (0.0,  0.0). 


Fig.  10(c)  Three-dimansional  display  at 
(0.0, 0.0)  point  of  view. 


Fig.  10(d)  Three-dimensional  display  at 
(0.0,0. 5)  point  of  view. 
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Fig.  11(a)  The  original  test  data. 


Fig.  11(b)  The  contour  map  with  JdB=1.0 
with  main  frequencies  around  ( 0 . 0 , 0 . 0 ) 
and  (0.06,0.125). 
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Fig.  11(c)  Three-dimensional  display  at 
(0.0, 0.0)  point  of  view. 


Fig.  11(d)  Three-dimensional  display  at 
(0.5, 0.5)  point  of  view. 


Fvm.  12(a)  The  original  test  data.  Fig.  12(b)  The  contour  map  with.idB=2. 

with  main  frequency’  around  (0.0 3,0.07). 


Fig.  12(c)  Three-dimensional  display  at  Fig.  12(d)  Three-dimensional  display  at 

(0.0, 0.0)  point  of  view.  (0.5, 0.5)  point  of  view. 
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Fin.  13(a)  The  test  data,  as  constructed  Fig.  -13 Ct>)  The  contour  map  withAdJ3=l. 

from  Fig.  11a  and  Fig.  12a).  with  main  frequency  around  (0.0, 0.0). 


Fig.  13(d)  Three-dimensional  display  at 
Fig.  13(c)  Three-dimensional  display  at  (0.5, 0.5)  roint  of  yiew. 

(0.0, 0.0)  point  of  view. 
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Fig.  14(c)  Tfiree-dirr.ensional  display  at 
(0.0, 0.0)  point  of  view. 


Fig.  14(d)  Three-dimensional  display  at 
(0.5, 0.5)  point  of  view. 
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Fir:.  15(a)  The  original  test  data.  Fig.  15(b)  The  contour  map  with  <ldB=-'3. 

with  main  frequency  around  (0.065,0.20) 


Fig .  16(c)  Three-dimensional  display  at 
(0.0, 0.0)  point  of  view. 


Fig.  16(d)  Three-dimensional  display  at 
(0.5, 0.5)  point  of  view. 
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The  two-dimensional  maximum  entropy  spectral  estimation  method  j 
proposed  by  Lim  and  Malik  is  studied  and  modified  for  spectral  i 
analysis  of  texture  images.  The  computational  efficiency  of  j 
the  algorithm  makes  it  feasible  for  minicomputer  implementation  ! 
to  extract  effective  spectral  domain  texture  features  for  i 
tex’'ture  discrimination.  Extensive  computer  results  are  [ 
presented.  i 
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